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We study finite size effects on ideal Bose gas thermodynamics (and where applicable Bose-Einstein 
condensation) in several systems with quadratic particle energy levels, namely, the Einstein universe, 
the one- dimensional box and the three-dimensional flat box (lDbox®R^). The Mellin- Barnes trans- 
form technique is used to obtain analytical approximations to thermodynamic quantities, providing 
an analytical description of the specific heat peak and the onset of condensation. 
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I. INTRODUCTION 



OO Since the experimental realization of Bose-Einstein condensation (BEC) for a dilute atomic gas in 1995, this field 

I ' has been the subject of intensive research, both theoretical and experimental. (See Q and for theoretical and 
experimental reviews respectively.) 

In particular, finite size effects have been dealt with in a number of works, both numerically and using a variety of 
analytical techniques. Rigorously speaking, the critical behaviour known as BEC may only happen when we take the 
limit N OO (N being the number of particles). For a finite system all thermodynamic quantities are analytical. As 
' the size of the system is taken to infinity in some well defined way (e.g., fixed temperature and density in the case of 
\^ . a homogeneous system), it is possible that the chemical potential (seen as a function of the temperature T when N is 
' given) will approach a non- analytical function at some critical temperature Tc. This will cause other thermodynamic 
,; quantities to also approach non-analytical behaviour, which can result in a phase transition. 

■ In this context, the harmonic oscillator trapping potential has received special attention |^-|^, motivated by its 
suitability to describe experimental traps and the simple form of the energy eigenvalues. The three-sphere (spatial 
section of the Einstein universe) has also received some attention. In this case the form of the energy eigenvalues is 
C \ still simple, but the fact that it involves squares of the quantum numbers describing the state makes the problem 
i of treating the thermodynamic sums more complicated. In fact, it is extremely complicated if the gas is considered 
^ ' to be relativistic. There is not to this date a study of the relativistic Bose gas in relying mostly on analytical 
O ■ expressions. Altaie [|j was the first to provide an analytical study of the non-relativistic case. Singh and Pathria |^ 
O studied the relativistic case obtaining Altaie's result in the non-relativistic limit. More recently, Trucks [|l^ repeated 
^ , the calculations of |^ for the minimally coupled case (whereas had done it for conformal coupling) and Altaie 
and Malkawi [|ll| for the spin-1 case (whereas |^ had done it for spin-0). Other authors have studied the problem of 
a Bose gas in an Einstein universe (and more generally, curved spaces) from different perspectives, not focusing on 
finite size effects ||T^[l5[|. 

Finite size effects for bosons trapped in a potential (including boxes) have also been studied analytically in a general 
setting [|l^. We will discuss this setting in more detail later. The basic idea behind the work of the Mellin-Barnes 
transform, is common to our work as well. There is nevertheless an important difference that makes the method of 
p6[ not applicable to the situations concerned here. 

All the investigations we have mentioned considered the particles to be non-interacting. Finite size effects have also 
been taken into account for bosons with self-interactions, but only numerically ]l7[ |. 

In this article, we make use of the Mellin-Barnes transform (used in similar contexts in the past in ||^, |l6| , |l8| , p^ ) to 
study the thermodynamics of an ideal non-relativistic spin-0 Bose gas in three systems that present some common 
features in the form of the energy eigenvalues: they are all quadratic. These systems are the three-sphere (non- 
relativistic spin-0 gas), the one-dimensional box (IDbox) and the IDbox^M^. (By this last example we mean a 
3-dimensional box, which is finite in one spatial direction and infinite in the other two.) 
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In section II, we use the Mellin-Barnes transform to obtain expressions for the number of particles and specific 
heat. In section III, we briefly discuss some matters inherent to this technique. In section IV, we expand the results 
of section II around the critical region to obtain the chemical potential as an explicit function of the temperature. 
From this, we obtain the specific heat as a function of the temperature only and consequently, we have an analytical 
description of the sharpening of the specific heat peak. In section V are some concluding remarks and a summary of 
the most important points in the article. 



II. NUMBER OF PARTICLES AND SPECIFIC HEAT 



Consider a system of non-self-interacting bosons with particle eigenstates of energy E„ ■ For a relatively high number 
of particles, and if we are not interested in fluctuations, which is the case, the thermodynamics of this system can be 
well described using grand-canonical statistics |^^|"p4|, greatly simplifying matters. The number of particles is then 

^ = ^[g/3(B„-p) _ ^ (1) 

n 

where P is the inverse of the temperature T, /x is the chemical potential and the sum is over all particle eigenstates. 
This expression gives us implicitly if we fix T and N. In this section, we shall aim at obtaining analytical expressions 
for TV and the specific heat, C, as functions of /i and (3. Rather than presenting the method in a general setting we 
prefer to develop it first for one specific case, the Einstein universe, and then generalize it to our other applications. 
We will use the natural units system (h = I and c = 1) and /cb = 1 throughout (/cb being Boltzmann's constant). 



A. The Einstein Universe 



In this case the energy eigenvalues are given by 

E,, = -\n^ + m^a^ + 6^-1]^ , n= 1,2,3,... (2) 
a 

with degeneracy c?„ = n^. a is the radius of the universe, m is the particle mass and ^ is the coupling constant. We 
will limit ourselves to conformal coupling, ^ = 1/6. As noted in [p] , this is not very restrictive since choosing a 
different coupling is equivalent to changing the mass from m to [rn^ + (6^ — l)/a^]^/^. For large a, this change is very 
small. In particular, in the large volume limit, all couplings will give the same results. In the non-relativistic limit we 
have 

Taking this limit requires ti^ ^ m?a?, ie, little occupation of high energy states ("high" when compared with to), 
corresponding to T <C to. For this to be true in the temperature range we are most interested in, the critical region, 
it suffices having p <C , where p is the particle density. 
We have then 



+ 00 



-i-oc 

n- y e ' , (4) 



n=l k=l 

where we have changed to the dimensionless variables, 



e = —2ma p . (5) 



2TOa^ ' 

(Note that e > —1 always and e —1 implies N — > -t-oo.) At this stage we invoke the Mellin-Barnes transform, 



1 

27ri 



c+ioo 

daT{a)x~°' , a: > 0, c> 0. (6) 
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Applying it to the exponential of (4) yields 



— / dar(a)2;-^fc-"^n2(n2+6)-", (7) 



TV 



where we have changed the order of summations and integration. For this we must have c > 3/2, so that the sums 
remain convergent. We now perform the integral by analytically continuing the sums to the whole complex plane 
(except eventual poles), closing the contour to the left in a large rectangle (corresponding to a high temperature, or 
small X, expansion) and applying the residue theorem. 



N ' 



-L f dar{a)x~^aa)h{a,e) . (8) 
27ri J 



7 is the contour mentioned above, C{a) is the Riemann zeta function and fi{a, e) is the analytical continuation of the 
sum X^ri^ n^'(n^ + e)~". The usefulness of defining fi{a,e) for general i becomes apparent when we consider other 
sums later in the article. These functions are treated in the appendix. We use the sign and not "=" as the part 
of the contour that was added to the integral is not necessarily vanishing. We will come to this later. 

In order to extract the poles of the integrand let us focus on r(Q;)/i(Q;, e). Consider for now a > 3/2. In this case 
/i(a, e) is just the summation we had before. By applying the F- function integral representation we have 

F(a)/i(a,e) = / ^-^e-*^ V Ti'e"*"'. (9) 

•^0 n=l 

We can now obtain an expansion for the sum in (9) by using the Mellin-Barnes transform on e~*" , yielding 

n2g-t«^ ^ — . / T{a)t-°'C{2a - 2) 



- f t-i , (10) 

where we solved the integral by closing the contour to the left and applying the residue theorem. For some properties 
of the F- and Riemann C- functions see Expression (9) now becomes 

F(a)/i(a,e) - ^ f dte-h-*' + ^ ^""^6"*% (11) 

where & is a positive number that can be arbitrarily small. The reason we divide the integral in two parts is that 
since the poles and their residues are the only pieces of information we require, we can concentrate on the first term 
of (11), which contains them. (The second integral is an analytic function of a.) We will denote the first term by A. 
Taylor expanding the exponential and performing the integral yields 

The analytical continuation to the complex plane is now trivial. The poles are at 

a=l~k, k = 0,l,2... (13) 



and the residues at these poles are 



Using this information in (0) results in 



resik)^^tf,^. (14) 



N - 4--^ ^C(? - kKex^ + x-'h{l,e) . (15) 



4 ^ fc! ^'2 

fc=0 
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The sum in this expression conies from the poles of r(a)/i(a, e), and the last term in comes from the pole of 
(^(a). In the high temperature range (small x) it is enough to consider only a few terms from the sum plus the other 
term. The value of /i(l,e) is given in (A14). We finally obtain the approximate result 



N . 



3th(7r\/i) 



(16) 



Instead of having a hyperbolic cotangent in the last term we could have it expressed in terms of a regular cotangent 
(see appendix). When we take the thermodynamic limit {V — > +00, constant density), we have x ^ Q and this result 
coincides with the continuum approximation, as would be expected. 

Altaic using a completely different technique (involving the Poisson summation formula), obtained a result 
almost identical to ours. The only difference, apart from the total disparity in the method, is that we have additional 
terms present. However, these new terms could have been found by extending the polylogarithm expansions in Altaie's 
method to higher order. 

We now turn our attention to the specific heat. The total internal energy of the system is given by 



1 



(17) 



where the sum is over particle energy eigenstates. In the case of non-relativistic particles in an Einstein universe, it is 

+00 



U ■ 



2ma 



1 



(18) 



The specific heat at constant radius of the universe is 

C 



dU 
dT 



N., 



5(3,1) 



5(2,1) 



5(1,1) 



where we define a class of sums 



\-oo +00 



5(i,j) = E"'^E^'^~'" 



(19) 



(20) 



n=l fc=l 



We used (^)jv = to obtain (^^) ■ (For details of how to get (19) see [||.) Applying the Mellin-Barnes 

V / N,a 

transform to the sums 1) results in 



Sit, f da r(a)x-"C(a " !)/.(«, e) , « = 1, 2, 3 , 



(21) 



as the analogue of (8). Studying T(a)fi{a, e) in the same manner as we did for the case z = 1, we find the poles given 
by 



1 



i~ k , fc = 0, 1,2, 



(22) 



The respective residues are 



res=-r 



2t + l\ , 
k\ " 



(23) 



Finally we have the expansion 
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+ 00 



fe=0 



k\ 



Using (A15), (A17) and (A18) for /i(2,e) yields 



5(1,1)-^ 
+ 

30F 



+ i-k] +fi{2,€)x- 



c ( ^ ) i-^r 



(24) 



-— = coth(7rv/e) -csch^(7r\/e) 

4ve ^ / 4 



5(2,1) 



5(3,1) 



_ 5 
a; 2 



+ 



- ^'^^ coth(7r-v/e) + ^csch^(7r-v/e) 
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a; 2 



57re2 



coth(7r-\/e) 



-csch^(7r-\/e) 



(25) 



where we consider xe small and truncate the sum in (24) at the third term. This can be easily extended if necessary 
to higher orders. 

We can now assess the accuracy of our approximation by comparing it with exact numerical results. We used (16) 
(with only the terms explicitly displayed) to determine e, from which we calculate the fraction of particles in the 
ground state given by Ng^/N = [e^(^+'^) — 1]~^/N and the specific heat from (19) and (25). Fig. 1 shows Ng,-/N for 
N = 10^ and N = 10^. Fig. 2 shows C/N for the same values of N. Had we also plotted the exact numerical results 
given by (19) and (20), the figures would look no different because the numerical plots would totally superimpose the 
analytical ones. For a better appreciation of the error involved in our analytical approximation see fig. 3, where we 
plot the error for A^gr. 

Some comments are in order here. First, in all figures the scaling factor for the temperature is the infinite volume 3D 
critical temperature, Tc, given by l/(2ma^Tc) = Xc = [{\/n(^{3/2))/{4N)]'^^^. Given the volume of the three-sphere, 
V = 27r^a^, we can put Tc in terms of the particle density, p, as Tc = {2tt /ra){p/ C,{2i/2))'^/^ . 

Secondly, we can clearly see from figs. 1 and 2 the approach to critical behaviour as N increases. In section IV we 
provide an analytical description for this phenomenon which explains the behaviour found here. 

Finally, we consider the errors seen in fig. 3. We can see that these are very small for a wide range of x around 
the critical region. The peak of ATYg^ = [iVg].(ana) — 7Vgr(num)]/Afgr(num) for TV — liT is a curious feature of our 
approximation. In a closer look, we can see this peak beginning to form already for N = 10^. It increases with N . 
This feature can be misleading. While it is unquestionable that AiVgr increases with N in the critical region, the error 
in (16) for a given e is actually diminishing. What happens is that as the thermodynamic limit is approached, we will 
have, for fixed T < Tc, e —1 and for fixed T > Tc, e +00, where Tc is the infinite volume critical temperature. 
This means that as N (and V) increases, e will be an increasingly steeper function of temperature at Tc. This causes 
e, given in (16) as an implicit function of N and x, to have an increasingly larger error in the critical region. The 
error is passed on to the ground state population. From fig. 3 we can also see that, excluding the already mentioned 
peak, the error increases as we get away from the centre (critical region). This is due to the series cut off at the 
third term in the relevant expressions. (|a;e| increases away from the centre.) By including more terms we can make 
the approximation work for an arbitrarily large range of x (or T). Similarly, neglecting more terms will cause the 
approximation to deteriorate faster as we get away from the centre. We do not plot the error in the specific heat, 
AC, because it does not show anything new. It behaves like ATVgr, except that the peak in the critical region is not 
manifest. This feature derives from C, as given in (19), being dominated by S{3, 1), which in turn is dominated by 
the term in a;~^/^ in (25). This term does not depend on e and so has a smooth behaviour. 



B. The One-Dimensional Box 



The energy eigenvalues for a non-relativistic spin-0 particle in a ID box are 



En 



2mL2 



n\ n = 1,2,3... 



(26) 
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where L is the length of the box. The levels are non-degenerate, ie, d„ = 1. After a change of variables, 
and using (1), we have 



(27) 



+ 00 



N 



J2 [e"^"'+^ 



n=l 



5(0,0) . 



Following the steps of the previous section yields 

1 



2^ / dar(a)x~"C(a)5(a,e) 



For T{a)g{a, e) the poles and respective residues are 

a = i - fc , A: = 0, 1,2, . 
2 



a 



-fc. A: = 0,1,2,. 



res = , , £ 

2 /c! 



res 



resulting in 



y +00 / -1 \fc 1 11 1 

fc=0 fe=l ^ ' 



(28) 
(29) 

(30) 
(31) 

(32) 



fc=0 

After truncating the sums and using (A3) we have 



N 



-X 2 



_l=coth(,rVi)-^ 



1 1 

4 " 24^ 



-ex 



For the specific heat we can easily derive 



C = x' 



5(2,1). 



5(1,1)2 



5(0,1) 

We already have analytical expressions for 5(2, 1) and 5(3, 1). 5(0, 1) is obtained in the same way, yielding 



(33) 



(34) 



5(0,1) 



/TT _ 1 
X 2 



An 



^ /5 
167r2^ I 2 



15 



1287r3 
1 



{exf 



— csch2(7rVi) + ^ coth(^V^) - — 



(35) 



where we use (AS). 

The behaviour of the specific heat and ground state population is uninteresting in the ID box case since all quantities 
are always smooth. It is well known ||2^,|2^ that for dimensions lower than three, there is no approach to BEC as 
N — * +CXD. For this reason, we do not include any plots here. The shapes of the curves can be easily extracted from 
the analytical expressions. In particular, in the large N limit we have 



N 



1 



^coth(vrVi)-^ 



(36) 



As to the error contained in the approximation, it is smaller in this case than in . So that the reader will 
have an idea of just how small it is, we will restrict ourselves to saying that for TV = 10^, AiVg^ = [iVgi.(ana) - 
A^gi (num)]/A'gi.(num) only becomes larger than 10~^ for x~^ — 10^. Furthermore, the larger the value of N , the 
smaller the error for fixed x or fixed x/xo, where xo is defined as e(a;o) = 0, in analogy to xo in 3D (see section IV). 
This happens because, unlike in the Einstein universe, e is bounded for any fixed x (or x/xq). 
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C. The ID Box (g) 



The lDbox(g)R^ can be defined in several ways. One of tliem is to consider a 3D rectangular box with sides of 
length Li, L2 and L3 and then keep Li fixed and take L2 = cL^ — > +00 with c e R+ being a constant. In a sense 
this is the most natural procedure and we shall adopt it. 

The energy eigenvalues of a non-relativistic particle in such a box are 



E„ 



TT 

2rn 



Using (1) and after some manipulation we have 

+00 +00 



-kf3 



2^2 



Hi 

L3 



1,2,3. 



+OC 77-^71^ 4-oc , 

e ^"-^-2 > e 



E 



Now, in the limit L9, 



k — l ni—l 712 — 1 713 — 1 

-00 the sums in 712 and can be replaced by integrals, yielding 



(37) 



(38) 



E 



-kl3- 



27rfc/3 



i-2,3, 



(39) 



as the leading terms. Inserting (39) in (38) and making the change of variables 

(3tt^ 2mLl 



2mLl ' 



we obtain 



N : 



^2^3 -1 



(40) 



(41) 



Since the volume is infinite when L2 and L3 are infinite, the number of particles is also infinite. It is then convenient 
to define a new quantity, 77, which will be useful later, as 



r;=-x-i5(0,-l) 



(42) 



•q is the number of particles in a cube with sides of length Li. We prefer using 77 instead of the particle density, as t] 
is the mathematical analogue of N in the previous cases. 

Similarly to the previous sections, we obtain an analytical approximation for S'(0, —1) by applying a Mellin-Barnes 
transform, yielding 



5(0,-1)^-^ / dar(a)a:-"C(a + l)5(a,e) . 



(43) 



Extracting the residues of the poles of the integrand is a little trickier in this case, because the pole at a = is double. 
This pole comes from T{a)Q{a + 1) since g{a^ e) is finite at a = 0. Furthermore, 



1 



r(a)C(a + l) = ^ + 0(a°) . 



(44) 



There is no term in a ^ (this is easily derived from Gamma and Riemann zeta function results in |p5| ). To get the 
residue we need the derivatives of x~°' and g{a, e) with respect to a. It will be 



d_ 

da 



3(a,e) - 5(0, e) Inx 



a=0 



Using (A9) we have for this pole 



1, , 2sinh(7rVi) 

res = - In a; — In r= 

2 



(45) 



(46) 
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The residues of the other poles are obtained in the manner of the previous sections, yielding 



1 



1 

4^^ 



1 2sinh(7rVi) 
+ - ma; — m ^= 

2 i/e 



(47) 



Like in the Einstein universe, this can be seen to approach the continuum approximation when x ^ 0. 
Now, for the specific heat, the analogue of expressions (19) and (34) here is 



C 



7rL2-^3 



x.(2.0)^2.(l,-l)^l.(0.-2)-x K"'°--;>;f"'°' 



(48) 



In the course of this article we already obtained analytical expressions for S{0, 0) (ID box number of particles), S{1, 0) 
{S^ number of particles) and 5(0, —1) (lDbox(g)M^ number of particles). We now need expressions for S{2, 0), 5(1,-1) 
and 5(0,-2). 

For 5(0, —2) the procedure is identical to the one used to get 5(0,-1). Now the double pole is at a = —1 and its 
residue is 



11, d 

res = —x€ xemx — — — 

2 2 da 



g{a,e) . 



(49) 



The expression for dg/da at a = —1, given by (All), is not as simple as we would have liked, as it includes an infinite 
sum. However, that turns out not to be limiting when it comes to studying analytically the critical region, as we will 
see. Using (All) we obtain 



5(0,-2)-^.-^ 



c 



1 1 , (,(3) , 2sinh(7rVi) 
+-ex — -exmx + -;r-^x + ea;ln ^ — x y 



27r2' 



^^(2 + 2fc)(2fc)! 



(50) 



for |e| < 1. i?2fe arc the Bernoulli numbers. 

5(2,0) is simpler. All the poles are simple and it does not present any additional complications in relation to the 
sums studied in the previous sections. The result is 



5(2,0) 



C 



1 . / I 



{exf 



3 

7re2 



coi\i{'n^/e)x ^ . 



(51) 



In the case of 5(1,-1), there are not any double poles either. Still, it deserves a little attention. Proceeding as 

before, we have 



5(1,-1)-^ j rfar(a)ar-"C(a + l)/i(a,e) . 



(52) 



The pole deriving from the zeta function is at a = 0. Despite the fact that r(a) has a pole at a = 0, the function 
T{a)fi{a, e) can be analytically continued there. The reason is that /i(0, e) = 0, cancelling the T{a) pole. The residue 
of the pole of r(a) at a = is 1. We have then. 



limr(a)A(a,e)= ^ 



/i(a,e) 



a=0 



C(3) +2g {2nf'^B,, 
^o(2 + 2A:)(2fc)! 



(53) 



The last equality is valid for |e| < 1 (see appendix). Finally we arrive at 



8 



^(1,-1)^^.- 



C(3) ^ (2.)^^i?.. 

27r2 ^ (2 + 2fc)(2fc)! ' ^ 

We can now obtain e from (42) and (47) and use it to calculate rii, which is the analog of Ng^. in the Einstein 
universe: 

7?l = -X^'^ ^ 



4 ^ fc 

fc=l 



4 L 



1 -e 



-x(l+e} 



(55) 



ryi is obtained from (42) by considering only the n = 1 part of the sum S'(0, — 1). It is not the ground state 
population. The fraction of particles in the ground state has no interest to us, since iVgr is finite (expression (38) with 
ni = n2 = = 1) and so, when L2 and L3 are infinite, Ng,-/N = 0. In contrast, rji/rj is the relevant quantity here. 
It has a similar behaviour to Ng,-/N in the Einstein universe. 

Making Li +00 in our system is a way of obtaining the infinite three dimensional space. And yet, we have 
Ngr/N — > for all temperatures when Li — > +00. This shows that the statement that in infinite 3D, Ngi-/N is 
finite below a certain critical temperature is incorrect. We must specify how we make our 3D space infinite. In 
the statement above, we should replace "infinite 3D" by, for example, "rectangular box of sides Li, L2 and L3 with 
Li,L2,L^ — > +00 while L1/L2 and £2/^3 are kept constant". One might rightly argue that the space we study here 
is not a box, since L2 and L3 are taken to infinity from the start. However, it is possible to show that we can have a 
3D box for which Ng,-/N — > for all temperatures when Li, L2, L3 — > +00 if the lengths of the sides obey particular 
relations between them. 

In fig. 1 we plot rii/r] for 77 = 10^ and ry ~ 10^. Only the analytical curves are plotted. Like in the Einstein 
universe, if the numerical curves were plotted, they would totally superimpose the analytical ones. In fig. 2 we plot 
C/N obtained from (48), (16), (33), (47), (50), (51) and (54). The point where the curves change the pattern in fig. 2 
is where e becomes larger than 1, rendering (50) and (54) not valid anymore. Both the analytical curves, obtained 
from the expressions mentioned above, and the exact numerical ones, obtained from (42), (48) and (20), are plotted. 
In the range where the analytical approximation applies (e < 1), they are totally superimposed, looking like one curve 
only. For e > 1 only the numerical curves are plotted. In both figs. 1 and 2, the curves are similar to those of the 
Einstein universe. In the large N limit they tend to be the same, specifically, the standard text book curves of the 
infinite box. This can be seen analytically by taking a; — *■ in our approximations. 

In fig. 3 we plot the error At^i — [771 (ana) — 77i(num)]/77i(num). All the error analysis done for the Einstein universe 
applies here also. 

III. A NOTE ON THE MELLIN-BARNES TRANSFORM 

As we mentioned already, when going from expression (7) to (8) we neglect the contribution from the added part 
of the contour to the integral. Although this contribution is vanishing for example for expression (6) ||2^ , it is not 
necessarily so for (7). It is at most very small, as we could see when comparing the analytical results to the exact 
numerical ones. However, the peak in fig. 3 clearly indicates that it is not vanishing. It is smaller and smaller as N 
gets larger at fixed density (and hence, x gets smaller), but if it was null the peak would not be there at all. 

It is possible, if this contribution is not negligible, to have a situation in which the Mellin-Barnes transform technique 
will give totally wrong results. In our investigations, we came across such an example. When applying the transform 
to (4) instead of applying it to the whole exponential, we could have split the exponential in two, e^*^^^" -i)g-kx(e+i) ^ 
take the ground state out |2^] , and apply the transform to the first part only, in the manner of ||l^ . This yields 

— / dar(a)a;-"Li„(e-^(^+i))/(«) , (56) 



where /(a) and the polylogarithm LiQ(a) are defined like 



+00 



/(a) = ^(n+l)2[n(ri + 2)]- 
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+ CXD 

fc=i 

Closing the contour to the left and proceeding as in the previous section, we obtain 

^-^^'^EFLii-^(«"^^^^'^)^'- (58) 
For e < 0, this expression is completely wrong. This can be seen by noting that 

+ 00 

E ^Lii-.(e-^(^+^^)2:'= = Lia (e^ ) . (59) 

Therefore, our expression is divergent for e < 0. Curiously, this is exactly the same as the continuum approximation 
(also called "bulk"), ie, directly replacing the sum in n in (4) by an integral. 

This way of applying the Mellin-Barnes transform is incorrect not only in the Einstein universe case. It is not 
difhcult to show that if the energy levels are of the form _E„ — _Eo = a(n + 6)* + c, a S R+, i G M"*", & S M \ Z~, c G K 
and the levels degeneracy, d„, is any polynomial (and this includes all the cases in this article), then the situation 
will be identical. In principle, there is nothing preventing this from generalizing to more cases still. Thus, one must 
be very careful about how to apply the Mellin-Barnes transform. 



IV. CRITICAL REGION EXPANSIONS 



Expressions (16) and (47) give us e implicitly and from here we have all thermodynamic quantities. However, it 
would be much more convenient to have e given explicitly. That is what we will aim for in this section, concentrating 
on the critical region. This is similar to what was done in in the context of a gas subject to a magnetic field. 
However, an important difference is that we do not perform our expansions around Xc- Instead, we use ccq (defined 
below), making the problem considerably easier to treat. This will still give us the critical behaviour since when 
N +00 we have xq — > Xc. We then use the results obtained on the specific heat. 



A. Einstein Universe 



Define xq as being the value of x at which e = 0. From (16), we have 



^-#cf^Vo--k-. (60) 



4 ' V2, 

Still in (16), change variables from x to 8 — [xjxi^)^!'^ — 1. The Taylor expansion of e as a function of 6 (for fixed A^) 



IS 



e = ai(5 + a2(5^ H , (61) 

where the OiS are, for the moment, unknown. There is no independent term since by construction e(0) — 0. We now 
expand coth(7r-ye) around -^e — and insert (60) and (61) in (16). Solving for each power of (5 finally yields 



i£c(i) 

For large iV, a;o ^ 1 and we have to highest order in x^^, 



ai^-x--- / . (62) 



ai~-%^Xo^. (63) 



The result for 02 is quite long and there is no point in displaying it in full. For xo ^ 1 it becomes 
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02 



20n 



(64) 



Expressions (61), (63) and (64) give us e explicitly in a neighbourhood of xq- 

We can now find an expansion for the specific heat by inserting (61) in (25), obtaining 



N 



15C(|) 

4c(i) 



+ 



81C(r _45C(|) 
407r 4C (I) 



729C(r 
1400V^ ° 



1.93- 1.38^- 



(65) 

The coefficients of each term are given only to highest order in Xq^ . The first term is the familiar specific heat 
maximum of the 3D infinite space. The other terms are dependent on which kind of 3D space we have, ie, they are 
specific to the Einstein universe. Prom here we immediately got the first and second derivatives of C with respect 
to 5 at (5 = and hence, with respect to T. The first derivative is a constant, but the second one goes to infinity as 
N +00 {xo 0). This is not surprising, it is the famous specific heat peak getting sharper. It does not mean 
that To is the same as the critical temperature T^. What we have is liniN^+ooTo = (which is already evident from 
(60)). However, if we define Ta as the temperature at which e = a 7^ 0, we will still have lim,{^^-\-ocTa = Tc because 
the function e(T) tends to become infinitely steep at T^. e = is, nevertheless, a privileged choice as it provides the 
simplest treatment for the problem. 

For completeness, we also present C/N at Tq expanded around xo =0 {N = +00): 



N^'^ 4C(i) 



+ 



15C(|) 27C(f) 



20FC(i)' 87r. 



~ 1.93 - 0.75x0' + 



(66) 



B. ID Box 



For the lDbox(8)M^ the process is exactly the same. In this case, xq is still the value of x at which e = and is 
given by (42) and (47) as 



0F^/3 
^ = — H2 



1 



Mxo) - -^Xq ^ ln(27r) . 



Again, we have that liniN- 
The result for e is 



>+<xXo 



9C(|) 

3 

7r2 



27Ci 



IOtt 



(67) 



(68) 



where S = \J xjxQ — 1. The coefficients arc only displayed to highest order in Xq ^ . 



For the specific heat, inserting (68) in (16), (33), (47), (50), (51) 



C_ 
N 



15C 



+ 



81C(i)' 45CI 



4c(i) 

1.93 + 3.02(5 



207r 



4C(|) 



6- 



and (54) we have 
243C(i)^ 



1400V7r 



1.75x0 



(69) 



Expressions (50) and (54) are valid only for |e| < 1. Since we are expanding around e = 0, this is not important, but 
it means that (69) is valid only for 5 > S' , where S' is a negative number satisfying e{S') = 1. Comparing (69) with 
(65) we see that the qualitative behaviour is very similar except that now the second term has a positive coefficient, 
meaning that Tq is to the right of the specific heat maximum, as opposed to the previous case. The second derivative 
still becomes infinitely negative as xq — > 0. Once more, this corresponds to the sharpening of the specific heat peak. 
The equivalent of (66) here is 



N 



(xo) 



15C(i) 15C(|) 



In 



Xo 



4C(|) 40FC(i)"° (2-) 



3c(i) 



27C(i) 

. 3 

477 2 



x^+- 



~ 1.93 - 0.42xo' Inxo - 2.35xo' + ■ 



(70) 
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V. CONCLUSION 



We have provided a general method to treat some thermodynamic problems by invoking the Mellin-Barnes trans- 
form. In principle, it can be used to treat any problem that presents the same kind of sums as long as the quantity 
equivalent to ex does not become very large in the range of interest. Although we have used it in systems with 
quadratic energy levels only, it should be possible to use it in other systems as well. This method is inspired by the 
one in with some important modifications. Moreover, we have shown that the method presented in p6| has limi- 
tations. In particular, it is not applicable to the spaces studied here and generally to systems with energy eigenvalues 
of some specific forms. Perhaps it would be interesting to look at this from the point of view of the density of states 
method which was shown in to be linked to the Mellin-Barnes technique. 

As already mentioned, the results obtained in section II for the number of particles in the Einstein universe could 
have been obtained using the Poisson summation formula as in |^ . It is likely that we could have used this technique 
to get all the results of section II. It would be interesting to establish a deeper connection between the two techniques. 

We have also provided an analytical description of the appearance of the sharp specific heat peak and hence, the 
onset of EEC. As far as we know, this has not been done before. It was a natural consequence of obtaining an 
analytical result for C, expanded around the critical region, in terms of the temperature only and not the chemical 
potential, which in turn, derived from having the chemical potential as an explicit function of the temperature, in a 
similar fashion to p8) . Our particular choice of temperature around which to perform the expansions simplified the 
calculations and therefore it was possible to obtain the results for the specific heat. 

As a side issue, from our treatment of the lDbox(g)M^ arose the conceptual problem of what one means by infinite 
three-dimensional space. We pointed out that the way this limit is taken must be specified before we can make certain 
statements about the system. In particular, the system treated is an example of the onset of EEC at a finite and 
not vanishing critical temperature while the fraction of particles in the ground state {Ngr/N) remains zero at all 
temperatures. Still, there will be a quantity that takes the standard role of N^^/N . In our case, it is the fraction of 
particles that has minimum frequency (ni = 1) in the Li direction. 
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APPENDIX: 



Consider the function defined, for i £ N*^, as 

-f oo 



/.(a,e) = + for 3fi(a) > z + ^ . (Al) 

When '^[a) < i -f 1/2, we will define /i(Q!, e) by analytic continuation as described below. In this appendix we want 
to obtain fi and dfi/da for some special values of i and a. 

Define g = Jq. The analytical continuation of this function can be achieved by noting that 

d 

— <7(a,e) = + l,e) . (A2) 

For 5R(a) > 1/2, this follows immediately from (Al). In order to have an analytical dg/de, this equality must extend 
to the complex plane. g{l, e) is given by 

+ 00 



^coth(.Vi)-^ (A3) 



or 
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5(1, e) 



(A4) 



Both expressions (A3) and (A4) are valid for e > and e < 0, but the first one is more convenient for e > and the 
second for e < 0. g{0, e) is obtained by integrating (A2) with a = 0, yielding 



We know that g{a, 0) = C{2a). Thus, 



<?(0,e)=ff(0,0) 



7(0,e) = C(0) = -- 



g{—l, e) is obtained by using (A2) with a = —1 in conjunction with (A6), resulting in 
We also need y(2, e). This comes directly from (A2). Putting = 1, we have 

TT^ TT 1 

5(2, e) = — csch2(7rVe) + coth(7r^/i) - — . 



(A5) 



(A6) 



(A7) 



(A8) 



Now, we want the derivatives dg/da at a = and a = —1. Inserting the Taylor series of g{a, e) and g(l + a, e) in 
(A2) and integrating, yields 



d_ 

da 



.9(a,e) = 2C'(0)- / g{l,e)de 



In 



2 sinh(7r-ye) 



(A9) 



where the prime in C'(0) means derivative and we used C'(0) = — (ln27r)/2 and (A3). Following the same procedure 
for dg/da at a = —1, we obtain 



d_ 

da 



g(a,e) = 2C'(-2) + ie- / In 
-1 ^ Jo 



2 sinh(7r-ye) 



de . 



(AlO) 



There is no simple form for the integral in (AlO). However, in the case of |e| < 1, it can be expressed in terms of an 
infinite sum (see eg pa]), which is more convenient than just leaving it as it is. After getting C'(~2) = — C(3)/(27r^) 
from expressions in p5[ , we finally have 



d_ 

da 



Q = -l 



C(3) , 2sinh(7rVi) 



k=0 



(2 + 2fc)(2fc)! 



(All) 



(valid for |e| < 1), where i?2fe are the Bernoulli numbers. 

The results for fi{a, e), i G N are now easy to obtain. These functions can be expressed in terms of g{a, e), as can 
be seen for 5R(q;) > i + 1/2 in (Al). For fi{a, e) we have 



/i(a, e) = g{a - 1, e) - eg{a, e) 



(A12) 



This equality must hold also for 3f?(a) < i + 1/2, so that fi(a, e) is analytical. The extension to i ^ 1 does not present 
any additional difficulties. The special values we will need are 



/i(0,e) = 



/i(l,e) 
/i(2,e) 



^ coth(7rVe) - Y*=sch^(^Ve) , 



/2(l,e) = ^coth(7rVi) 



(A13) 
(A14) 

(A15) 
(A16) 
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d_ 

da 



/3(2,e) 
/i(a,e) 



a=0 



coth(7r-ye) 



-csch^(7rye) , 



^ coth(7r^/e) ^csch^(7r-v/e) , 



C(3) 
27r2 



E 

fc=0 



(2 + 2fc)(2fc)! 



(A17) 
(A18) 
(A19) 



Expression (A19) is valid only for |e| < 1. 



Y. Castin, in "Coherent Matter Waves", Lecture Notes of Les Houches 1999 Summer School, Session LXXII, edited by 
R. Kaiser, C. Westbrook and F. David (EDP Sciences and Springer- Verlag, 2001); F. Dalfovo, S. Giorgini, L. P. Pitaevskii 
and S. Stringari Rev. Mod. Phys. 71, 463 (1999). 

W. Ketterle and E. Cornell, in "Coherent Matter Waves", Lecture Notes of Les Houches 1999 Summer School, Session 

LXXII, edited by R. Kaiser, C. Westbrook and F. David (EDP Sciences and Springer- Verlag, 2001). 

S. Grossmann and M. Holthaus, Phys. Lett. A 208, 188 (1995). 

W. Ketterle and N. J. van Druten, Phys. Rev. A 54, 656 (1996). 

K. Kirsten and D. J. Toms, Phys. Lett. A 222, 148 (1996). 

K. Kirsten and D. J. Toms, Phys. Rev. A 54, 4188 (1996). 

H. Haugerud, T. Haugset and F. Ravndal, Phys. Lett. A 225, 18 (1997); T. Haugset, H. Haugerud and J. O. Andersen, 
Phys. Rev. A 55, 2922 (1997). 

M. B. Altaie, J. Phys. A: Math. Gen. 11, 1603 (1978). 

S. Singh and R. K. Pathria, J. Phys. A: Math. Gen. 17, 2983 (1984). 

M. Trucks, Phys. Lett. B 445, 117 (1998). 

M. B. Altaie and E. Malkawi, J. Phys. A: Math. Gen. 33, 7093 (2000). 
L. Parker and Y. Zhang, Phys. Rev. D 44, 2421 (1991). 

D. J. Toms, Phys. Rev. Lett. 69, 1152 (1992); D. J. Toms, Phys. Rev. D 47, 2483 (1993). 
K. Kirsten and D. J. Toms, Phys. Rev. D 51, 6886 (1995). 
J. D. Smith and D. J. Toms, Phys. Rev. D 53, 5771 (1996). 

K. Kirsten and D. J. Toms, Phys. Lett. A 243, 137 (1998); K. Kirsten and D. J. Toms, Phys. Rev. E 59, 158 (1999). 

F. Brosens, J. T. Devreese and L. F. Lemmens, Solid State Comm. 100, 123 (1996). 

G. B. Standen and D. J. Toms, Phys. Rev. E 60, 5275 (1999). 

M. Holthaus, E. Kalinowski and K. Kirsten, Annals of Physics 270, 198 (1998). 

H. D. Politzer, Phys. Rev. A 54, 5048 (1996). 

M. Gajda and K. Rzazewski, Phys. Rev. Lett. 78, 2686 (1997). 
M. Wilkens and C. Weiss, J. Mod. Opti cs 44, 1801 (1997) 
S. Grossmann and M. Holthaus, e-print cond-mat /970904^ 



P. Borrmann, J. Harting, O. Mulken and E. R. Hilf, Phys. Rev. A 60, 1519 (1999). 
L S. Gradshteyn and L M. Ryzhik, Table of Integrals Series and Products (Academic Press, 1994). 
R. M. May, Phys. Rev. 115, 254 (1959). 
K. Kirsten and D. J. Toms, Phys. Lett. B 368, 119 (1996). 
In fact, this is how you prove (6) (see Q). 

The reason we do the split like this and not e~*^" e'''^" and take the ground state out is to allow the separation of the 
sums in k and n, with each of them remaining convergent. 



14 




0.0 



1 2 



FIG. 1. Analytical Ngr/N in the Einstein universe for — 10^ (continuous line) and A'^ = 10^ (dotted line) and rii/r] in the 
lDbox®R^ for rj = 10'^ (dashed line) and rj — 10^ (dashed/dotted line) as functions of the rescaled temperature 
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FIG. 2. The specific lieat per particle in tlie Einstein universe for A*' — 10^ (analytical; dashed line) and A'' = 10^ (analytical; 
dashed/dotted line) and in the lDbox®R^ for rj — 10^ (continuous -analytical- and then dashed -numerical- line) and 77 = 10^ 
(dashed -analytical- and then dotted -numerical- line) - see text for details 
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FIG. 3. AA'gr in the Einstein universe for A'^ = 10^ (continuous line) and A*' = 10^ (dotted line) and Aj^i in the lDbox®R^ 
for r] — 10^ (dashed line) and rj = 10^ (dashed/dotted line) - see text for the definitions of AA'gr and Arji 
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